clear; clc; close all;

% 参数设置
a   = -1;       % 忆阻器参数 a
b0  = 2;        % 忆阻器参数 b 的某个默认值
q0  = 0;        % 初始内部状态
dt  = 1;        % 步长 \Delta t
N   = 100;       % 总迭代步数 (可自行调整)
nVals = 0:N-1;  % 用于计算 i_n, v_n 的离散时刻


figure('Position',[100,100,1200,400]);
sgtitle('离散忆阻器在的轨迹');

% (a) 改变 \omega = 0.1, 0.2, 0.3
subplot(1,3,1);
hold on; box on; grid on;

% 固定 F=0.2, a=-1, b=2, q0=0
F_a  = 0.2;
b_a  = 2;
omega_list_a = [0.1, 0.2, 0.3];
cmap = lines(numel(omega_list_a));  % 生成颜色

for k = 1:numel(omega_list_a)
    w = omega_list_a(k);
    
    % 初始化
    q = q0;             % q_n
    i_data = zeros(1,N);
    v_data = zeros(1,N);

    for n = 1:N
        in   = F_a * sin(w*(n-1));      % i_n
        vn   = (a + b_a*q^2) * in;      % v_n
        i_data(n) = in;
        v_data(n) = vn;
        % 更新内部状态
        q = q + dt*in;  % q_{n+1}
    end

    % 绘制 (i_n, v_n) 轨迹
    plot(i_data, v_data, 'Color', cmap(k,:), 'LineWidth',1.5, ...
         'DisplayName', ['\omega=', num2str(w)]);
end

xlabel('i_n'); ylabel('v_n');
title('(a) F=0.2, q_0=0, a=-1, b=2');
legend('Location','best');
hold off;


% 改变 F = 0.1, 0.2, 0.3 ==============
subplot(1,3,2);
hold on; box on; grid on;

% 固定 \omega=0.1, a=-1, b=2, q0=0
w_b    = 0.1;
b_b    = 2;
F_list = [0.1, 0.2, 0.3];
cmap = lines(numel(F_list));

for k = 1:numel(F_list)
    F_val = F_list(k);
    
    % 初始化
    q = q0;
    i_data = zeros(1,N);
    v_data = zeros(1,N);

    for n = 1:N
        in   = F_val * sin(w_b*(n-1));
        vn   = (a + b_b*q^2) * in;
        i_data(n) = in;
        v_data(n) = vn;
        % 更新内部状态
        q = q + dt*in;
    end

    plot(i_data, v_data, 'Color', cmap(k,:), 'LineWidth',1.5, ...
         'DisplayName', ['F=', num2str(F_val)]);
end

xlabel('i_n'); ylabel('v_n');
title('(b) \omega=0.1, q_0=0, a=-1, b=2');
legend('Location','best');
hold off;


% 改变 b = 1, 2, 3
subplot(1,3,3);
hold on; box on; grid on;

% 固定 \omega=0.1, F=0.2, a=-1, q0=0
w_c  = 0.1;
F_c  = 0.2;
b_list = [1, 2, 3];
cmap = lines(numel(b_list));

for k = 1:numel(b_list)
    b_val = b_list(k);
    
    % 初始化
    q = q0;
    i_data = zeros(1,N);
    v_data = zeros(1,N);

    for n = 1:N
        in   = F_c * sin(w_c*(n-1));
        vn   = (a + b_val*q^2) * in;
        i_data(n) = in;
        v_data(n) = vn;
        % 更新内部状态
        q = q + dt*in;
    end

    plot(i_data, v_data, 'Color', cmap(k,:), 'LineWidth',1.5, ...
         'DisplayName', ['b=', num2str(b_val)]);
end

xlabel('i_n'); ylabel('v_n');
title('(c) \omega=0.1, F=0.2, q_0=0, a=-1');
legend('Location','best');
hold off;
